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A  significant  difference  between  the  behavior  in  tension  versus  compression  is  obtained  at  the  polycrystal 
level  if  either  twinning  or  non-Schmid  effects  are  contributors  to  the  plastic  deformation  at  the  single 
crystal  level.  Examples  of  materials  that  exhibit  tension-compression  asymmetry  include  hexagonal 
close-packed  (HCP)  polycrystals  and  intermetallics  (e.g.,  molybdenum  compounds).  Despite  recent  pro¬ 
gress  in  modeling  their  yield  behavior  in  the  absence  of  voids,  the  description  of  coupling  between  plas¬ 
ticity  and  damage  by  void  growth  in  these  materials  remains  a  challenge. 

This  paper  is  devoted  to  the  development  of  a  macroscopic  anisotropic  yield  criterion  for  a  porous 
material  when  the  matrix  material  is  incompressible,  anisotropic  and  displays  tension-compression 
asymmetry.  The  analytical  yield  criterion  is  obtained  based  on  micromechanical  considerations  and 
non-linear  homogenization.  The  matrix  plastic  behavior  is  described  by  the  Cazacu  et  al.  (2006)  aniso¬ 
tropic  yield  criterion  that  is  pressure-insensitive  and  accounts  for  strength-differential  effects.  Compar¬ 
ison  between  finite  element  cell  calculations  and  theory  show  the  predictive  capabilities  of  the  developed 
anisotropic  model  in  terms  of  modeling  the  combined  effects  of  anisotropy,  tension-compression  asym¬ 
metry  of  the  matrix  and  voids  on  the  overall  yielding  of  the  porous  aggregate.  It  is  shown  that  if  the 
matrix  material  does  not  display  tension-compression  asymmetry,  the  developed  criterion  reduces  to 
that  of  Benzerga  and  Besson  (2001 ).  If  the  matrix  is  isotropic,  it  reduces  to  the  isotropic  criterion  devel¬ 
oped  in  Cazacu  and  Stewart  (2009). 

©  2010  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

Ductile  failure  in  metals  occurs  due  to  the  nucleation,  growth 
and  coalescence  of  voids  (McClintock,  1968;  Rousselier,  1987). 
Voids  are  nucleated  in  metals  mainly  by  decohesion  at  the  parti¬ 
cle-matrix  interfaces  or  by  micro-cracking  of  second-phase  parti¬ 
cles  (see,  for  example,  Tvergaard,  1981).  Additionally,  voids  can 
nucleate  in  single  crystals  that  contain  neither  pre-existing  voids 
nor  inclusions  (see,  for  example,  Cuitino  and  Ortiz,  1996;  Lubarda 
et  al.,  2004  and  the  recent  studies  on  cylindrical  void  growth  in  ri- 
gid-ideally  plastic  single  crystals  of  Kysar  et  al.,  2005,  2006;  Kysar 
and  Gan,  2007).  Thus,  the  ability  to  accurately  describe  the  evolu¬ 
tion  of  voids  in  a  ductile  metal  is  crucial  to  being  able  to  accurately 
predict  its  failure. 

Gurson  (1977)  developed  widely  used  macroscopic  yield  crite¬ 
ria  for  porous  metals  containing  either  spherical  or  cylindrical 
voids  and  with  the  matrix  obeying  the  von  Mises  isotropic  yield 
condition.  The  success  of  Gurson’s  (1977)  criterion  lies  in  the  fact 
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that  it  was  deduced  based  on  micromechanical  considerations. 
Several  modifications  of  Gurson’s  (1977)  criterion  were  proposed 
(see  Tvergaard,  1981;  Tvergaard  and  Needleman,  1984;  Koplik 
and  Needleman,  1988)  based  on  finite  element  calculations  to  ac¬ 
count  for  void  interaction  and  coalescence.  Gologanu  et  al.  (1993, 
1994,  1997)  generalized  Gurson’s  (1977)  analysis  by  considering 
a  spheroidal  volume  containing  a  confocal  spheroidal  cavity.  In 
Garajeu  (1995)  and  Garajeu  et  al.  (2000)  the  overall  response  of  a 
porous  metal  with  a  viscous  matrix  was  investigated.  Gurson’s 
analysis  has  also  been  extended  to  the  case  when  the  matrix  mate¬ 
rial  is  compressible  and  obeys  a  Drucker-Prager  yield  criterion  (see 
Jeong  and  Pan,  1995;  Guo  et  al.,  2008). 

Most  metallic  alloys  display  plastic  anisotropy  as  a  result  of 
forming  processes.  Recent  studies  have  been  devoted  to  the  exper¬ 
imental  characterization  of  the  anisotropy  of  fracture  in  different 
alloys  (see,  for  example,  Benzerga  et  al.,  2004a,  for  an  overview 
of  experimental  evidence  for  anisotropic  ductile  fracture  in  steels). 
Liao  et  al.  (1997)  extended  Gurson’s  (1977)  cylindrical  criterion  to 
account  for  transverse  isotropy  by  using  Hill’s  (1948)  yield  crite¬ 
rion  for  the  matrix  material;  however,  these  authors  assumed  that 
the  anisotropy  in  the  plane  of  the  sheet  is  weak  and  can  be  de¬ 
scribed  by  a  single  anisotropy  coefficient.  Benzerga  and  Besson 
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(2001)  extended  Gurson’s  (1977)  criterion  for  fully  orthotropic 
metals.  Upper-bound  yield  surfaces  for  both  a  hollow  sphere  and 
cylinder  were  deduced  by  assuming  that  the  matrix  material  can 
be  described  using  Hill’s  (1948)  criterion  with  six  non-zero 
anisotropy  coefficients  for  three-dimensional  stress  conditions. 
Assuming  triaxial  loading  conditions  aligned  with  the  material 
symmetry  axes,  closed-form  approximate  yield  criteria  were 
obtained  for  both  spherical  and  cylindrical  voids. 

Numerical  studies  have  been  conducted  to  verify  and  validate 
these  anisotropic  Gurson-like  criteria  (e.g.,  Chien  et  al.,  2001; 
Wang  and  Pan,  2004).  Generally,  finite-element  analyses  of  a  cube 
containing  a  spherical  void,  subjected  to  plane-stress  conditions, 
were  performed  for  various  initial  porosities  and  different  values 
of  the  anisotropy  parameters.  More  recently,  these  dilatational 
anisotropic  plastic  models  have  been  used  to  predict  forming  limits 
for  anisotropic  sheets  containing  voids  (see,  for  example,  Brunet 
et  al.,  2004). 

Note  that  the  anisotropic  models  cited  previously  can  describe 
the  response  of  porous  media  only  for  particular  void  shapes 
(spherical  or  cylindrical)  and  axisymmetric  loading  conditions.  It 
should  be  expected  that  void  shape  has  a  significant  influence  on 
the  behavior  of  porous  anisotropic  metals,  yet  there  have  been  rel¬ 
atively  few  studies  on  the  combined  effects  of  void  shape  and  tex¬ 
ture  (see  for  example  the  experimental  study  of  Benzerga  et  al. 
(2004a)).  Benzerga  et  al.  (2004b)  proposed  a  yield  criterion  which 
combines  properties  from  both  Gologanu  et  al.’s  (1993)  criterion 
and  Benzerga  and  Besson’s  (2001)  criterion  to  account  for  both 
void  shape  (elliptical)  and  orthotropy,  respectively.  Recently,  Mon- 
chiet  et  al.  (2008)  used  a  limit  analysis  approach  to  develop  an  ana¬ 
lytical  yield  criterion  for  a  matrix  material  obeying  Hill’s  (1948) 
anisotropic  yield  criterion  and  containing  elliptical  voids.  The  limit 
analysis  was  performed  for  arbitrary  deformation  of  the  represen¬ 
tative  volume  element.  Furthermore,  it  was  shown  that  for  the  case 
when  the  matrix  is  isotropic,  Monchiet  et  al.’s  (2008)  criterion  pro¬ 
vides  a  rigorous  generalization  of  the  Gologanu  et  al.  (1993)  model. 
In  Keralavarma  and  Benzerga  (2010)  further  investigation  into  the 
evolution  of  the  voids’  orientation  was  provided. 

All  the  studies  cited  so  far  assume  that  the  yield  strengths  in 
tension  and  compression  are  equal  in  the  void-free,  or  matrix, 
material.  However,  in  the  absence  of  voids,  some  materials  with 
cubic  crystal  structure  (see  Benzerga  et  al.,  2004a,  for  experimental 
data  on  high  strength  steels)  are  pressure-insensitive  but  exhibit 
tension-compression  asymmetry  in  yielding.  This  shear  related 
strength-differential  (S-D)  effect  can  result  from  single  crystal 
plastic  deformation  due  either  to  twinning  or  to  slip  that  does 
not  obey  the  well-known  Schmid  law  (see,  for  example,  Vitek 
et  al.,  2004;  Hosford  and  Allen,  1973).  Cazacu  and  Stewart  (2009) 
recently  extended  Gurson’s  (1977)  analysis  to  the  case  when  the 
matrix  material  is  incompressible  but  displays  tension-compres¬ 
sion  asymmetry  and  developed  an  isotropic  plastic  potential  for  a 
porous  aggregate  that  is  sensitive  to  the  third  invariant  of  the 
stress  deviator. 

Metals  with  hexagonal  crystal  structure  exhibit  both  ten¬ 
sion-compression  asymmetry  and  pronounced  anisotropy.  Experi¬ 
mental  studies  have  shown  that  failure  occurs  by  void  growth  and 
coalescence  (see  Huez  et  al.,  1998).  A  fundamental  issue  that  arises 
is  how  the  texture  and  tension-compression  asymmetry  of  the  ma¬ 
trix  influences  the  void  growth  stage  of  the  failure  process.  A  key 
ingredient  in  the  development  of  a  ductile  failure  model  is  an  effec¬ 
tive  yield  criterion.  The  aim  of  this  paper  is  to  develop  a  closed- 
form  macroscopic  yield  criterion  for  anisotropic  porous  aggregates 
containing  spherical  voids  using  a  micromechanical  approach. 

The  structure  of  this  article  is  as  follows.  In  Section  2,  the 
homogenization  approach  due  to  Hill  (1967)  and  Mandel  (1972) 
that  is  used  in  developing  the  macroscopic  plastic  potential  for 
the  void-matrix  aggregate  is  presented.  Next,  the  anisotropic 


Cazacu  et  al.  (2006)  yield  criterion  being  used  to  model  the  rigid- 
plastic  behavior  of  the  matrix  material  in  the  void-matrix  aggre¬ 
gate  is  recalled  (see  Section  3).  Expressions  for  the  macroscopic 
stresses  and  overall  plastic  dissipation  are  derived  in  Section  4. 
The  main  result  of  this  paper,  the  expression  of  the  analytical 
anisotropic  plastic  potential,  is  presented  in  Section  5.  Finally, 
the  accuracy  of  the  developed  criterion  for  plastic  anisotropic  med¬ 
ia  displaying  strength-differential  effects  is  assessed  through  com¬ 
parison  with  numerical  results  obtained  using  finite  element  cell 
calculations  (see  Section  6). 


2.  Kinematic  homogenization  approach  of  Hill  and  Mandel 


Consider  a  representative  volume  element  Q ,  composed  of  a 
homogeneous  rigid-plastic  matrix  and  a  traction-free  void.  The 
matrix  material  is  described  by  a  convex  yield  function  cp{a)  in 
the  stress  space  and  an  associated  flow  rule 


a) 


where  a  is  the  Cauchy  stress  tensor,  A  ^  0  denotes  the  plastic  mul¬ 
tiplier  rate  and  d  =  (l/2)(  Vv  +  VvT)  denotes  the  rate  of  deformation 
tensor  with  v  being  the  velocity  field.  The  yield  surface  is  defined  as 
cp(t j)  =  0.  Let  C  denote  the  convex  domain  delimited  by  the  yield 
surface  such  that 


C=  {<%(*)<()}.  (2) 

The  plastic  dissipation  potential  of  the  matrix  is  defined  as 
w(d)  =  sup(<7  :  d),  (3) 

creC 


where  denotes  the  tensor  double  contraction.  Uniform  rate  of 
deformation  boundary  conditions  are  assumed  on  the  boundary  of 
the  RVE,  0£2,  such  that 

v  =  D  x  for  any  x  e  dQ  (4) 

with  D,  the  macroscopic  rate  of  deformation  tensor,  being  constant 
and  x  being  the  Cartesian  position  vector.  For  the  boundary  condi¬ 
tions  given  by  Eq.  (4),  the  Hill  (1967),  Mandel  (1972)  lemma  ap¬ 
plies;  hence, 

<ff:d>fl  =  E:D,  (5) 


where  ( )  denotes  the  average  value  over  the  representative  volume 
Q ,  and  £  =  {a) a-  Furthermore,  there  exists  a  macroscopic  plastic 
dissipation  potential  W(D)  such  that 


dW(D) 

~dD~ 


(6) 


with 


W(  D) 


inf  (w(d)) 

deK(D) ' 


Q’ 


(7) 


where  K(D)  is  the  set  of  incompressible  velocity  fields  satisfying  Eq. 
(4)  (for  more  details  see  Gologanu  et  al.,  1997;  Leblond,  2003).  This 
result  will  be  further  used  to  derive  the  plastic  potential  of  the  void- 
matrix  aggregate  in  the  case  of  a  random  distribution  of  spherical 
voids.  To  account  simultaneously  for  anisotropy  and  tension-com¬ 
pression  asymmetry  associated  with  directional  shear  mechanisms 
at  the  single  crystal  level,  the  matrix  material  is  assumed  to  obey 
the  Cazacu  et  al.  (2006)  anisotropic  yield  criterion.  This  yield  crite¬ 
rion  and  its  dual  in  the  strain  rate  space  will  be  briefly  presented 
next. 
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3.  Anisotropic  yield  criterion  for  the  matrix  material 

To  account  for  both  strength  differential  effects  and  anisotropy 
in  pressure-insensitive  hexagonal  metals,  Cazacu  et  al.  (2006)  pro¬ 
posed  a  3-D  orthotropic  yield  criterion  (denoted  in  the  following  as 
CPB06).  This  criterion  is  an  extension  to  orthotropy  of  the  follow¬ 
ing  isotropic  yield  function: 

G(s1,s2,s3; /<, a)  =  (|S!|  -  ks,)a  +  (|s2|  -  ks2)a  +  (|s3|  -  ks3)a  (8) 

where  s3,  s2  and  s3  are  the  principal  values  of  the  stress  deviator.  The 
material  parameter  k  captures  strength  differential  effects  while  a  is 
the  degree  of  homogeneity.  Starting  from  the  isotropic  function  gi¬ 
ven  by  Eq.  (8),  anisotropy  is  then  introduced  through  a  linear  trans¬ 
formation  operating  on  the  Cauchy  stress  deviator  a'.  The  general 
form  of  the  CPB06  anisotropic  yield  function  is  as  follows: 

F(&u&2, &3;k, a)  =  (|ff!  |  -  /cff,)0  +  (|ff2|  -  fcff2)°  +  (|ff3|  -  ka3)a, 

(9) 

where  <7i,<72  and  cr3  are  the  principal  components  of  the  trans¬ 
formed  stress  tensor 


with  material  parameters  o\  and  m.  As  previously  stated,  the 
parameter  o\  is  the  uniaxial  tensile  yield  strength  along  an  axis  of 
orthotropy  while  the  constant  m  is  defined  such  that  Eq.  (14)  is 
identically  satisfied  for  uniaxial  tensile  loading  along  this  orthotro¬ 
py  axis. 

Let  (ei,  e2,  e3)  be  the  reference  frame  associated  with  orthotro¬ 
py.  In  the  case  of  a  plate,  ei,  e2  and  e3  represent  the  rolling,  trans¬ 
verse  and  through-thickness  directions,  respectively.  Relative  to 
the  orthotropy  axes,  the  fourth-order  tensor  L  (see  Eq.  (10))  is  rep¬ 
resented  in  Voigt  notation  by 
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The  constant  m  is  expressible  in  terms  of  the  anisotropy  coeffi¬ 
cients  Ltj  and  the  strength  differential  parameter  k  as 


<7  =  L:<7'.  (10) 

In  Eq.  (10),  a'  is  the  Cauchy  stress  deviator  and  L  is  a  fourth-order 
symmetric  tensor  invariant  with  respect  to  the  orthotropy  group. 
Thus,  the  yield  condition  can  be  written  as 

ffe  —  ffi>  (11) 

where  oe  is  the  effective  stress  associated  to  the  yield  function  of  Eq. 
(9)  and  o\  is  the  uniaxial  tensile  yield  strength  along  an  axis  of 
orthotropy  (e.g.,  the  rolling  direction).  In  other  words, 

ffe  =  ce  [(|ffi|  -  ka i)“  +  (|ff2|  -  ka2)a  +  (|ff3|  -  ka3)a]v“,  (12) 

where  ce  is  a  constant  defined  such  that  oe  reduces  to  the  tensile 
yield  stress  along  an  axis  of  orthotropy. 

In  Cazacu  et  al.  (2006),  the  physical  significance  of  the  coeffi¬ 
cients  involved  in  the  CPB06  yield  function  given  by  Eq.  (8)  was  pre¬ 
sented  and  an  identification  procedure  based  on  the  results  of  tensile 
and  compressive  tests  was  outlined.  It  has  been  shown  that  the 
orthotropic  CPB06  yield  function  given  by  Eq.  (9)  describes  with 
accuracy  the  yield  loci  of  various  hexagonal  materials  (see,  for  exam¬ 
ple,  Cazacu  et  al.,  2006;  Plunkett  et  al.,  2008;  Khan  et  al.,  2007).  For 
the  magnesium  alloys  and  zirconium  materials,  a  value  of  the  homo¬ 
geneity  parameter  a  =  2  (see  Eq.  (9))  approximates  the  plastic  behav¬ 
ior  best.  Thus,  in  this  paper  a  =  2  will  be  used  in  Eq.  (9)  to  describe  the 
yield  behavior  of  the  matrix.  It  is  worth  noting  that  if  the  material 
does  not  display  tension-compression  asymmetry,  then  the 
strength-differential  parameter  k  is  automatically  equal  to  zero.  If 
k  =  0,a  =  2  and  L  is  constrained  to  be  deviatoric,  the  CPB06  yield  cri¬ 
terion  reduces  to  that  of  Hill  (1948). 


3.1.  Local  stress  potential 


For  a  =  2,  the  effective  stress  of  Eq.  (12)  becomes 


Ge  =  m 


\ 


-  kerf. 


(13) 


The  constant  m  is  defined  such  that  oe  reduces  to  the  tensile  yield 
stress  along  a  specified  axis  of  orthotropy.  Due  to  the  assumption 
of  associated  plastic  flow,  the  stress  potential  for  the  matrix  can 
be  written  as 


(16) 


with 

,  2  1.  1. 

(P\  =  ^  Mi  -  -Li2  -  -Li3. 

@2  =  12  3  ^22  3  ^23-  (17) 

@3  =  13  —  ^ f 23  —  ^  ^33- 

In  the  same  spirit  as  Hill  (1948),  it  is  assumed  here  that  the  sum 
of  the  L-components  on  the  first  row  is  equal  to  the  sum  of  the  L- 
components  on  the  second  row  and  to  the  sum  of  the  L-compo¬ 
nents  on  the  third  row.  Specifically, 

1-11  +  k\2  +  ^13  =1, 

L12  +  L22  +  £23  =  1,  (18) 

Li3  +  L23  +  L33  =  1 . 


The  constant  on  the  right  hand  side  of  Eq.  (18)  has  been  chosen  to 
be  unity  such  that  for  isotropic  conditions,  the  tensor  L  reduces  to 
the  fourth-order  identity  tensor,  1.  Note  that  the  additional  con¬ 
straints  of  Eq.  (18)  ensure  that  the  transformed  stress  tensor  d  is 
deviatoric.  As  an  example,  in  Fig.  1  is  shown  the  representation  in 
the  biaxial  plane  (an,  a33)  of  the  yield  surface  given  by  Eq.  (14)  cor¬ 
responding  to  several  materials  displaying  tension-compression 
asymmetry  characterized  by  the  same  value  of  the  strength-differ¬ 
ential  parameter  k  =  0.3098.  One  of  these  materials  is  isotropic 
while  the  other  two  materials  are  transversely  isotropic  with  (ei, 
e2)  being  the  plane  of  symmetry  (one  material  is  stronger  in  the 
plane  of  symmetry  while  the  other  one  is  stronger  in  the  direction 
normal  to  the  symmetry  plane,  e3). 


3.2.  Local  plastic  dissipation 

A  key  step  in  obtaining  the  overall  plastic  dissipation  W(D)  of 
the  void-matrix  aggregate  is  the  calculation  of  the  local  plastic  dis¬ 
sipation  w(d)  (see  Eq.  (3))  associated  with  the  matrix  yield  crite¬ 
rion.  Since  the  anisotropic  criterion  that  describes  yielding  in  the 
matrix  is  a  homogeneous  function  of  first  order  in  stresses  (see 
Eq.  (14)),  the  local  plastic  dissipation  becomes 


(p(d-m,G\)  =  m 


^(|ffi|-kffl)2-ff[  =  0 


(14) 


w(d)=A<r[.  (19) 

In  other  words,  A  is  the  dual  of  the  anisotropic  stress  potential  cp  gi¬ 
ven  by  Eq.  (14). 
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Fig.  1.  Plane  stress  yield  loci  for  void-free  materials  according  to  the  CPB06  yield 
criterion  given  by  Eq.  (14).  The  solid  line  denotes  the  isotropic  material  with  the 
other  curves  represent  the  two  transversely  isotropic  materials.  All  these  materials 
display  tension-compression  asymmetry  with  k  =  0.3098  (i.e.,  the  tensile  yield 
strengths  are  greater  than  the  compressive  yield  strengths). 


In  the  (ei,  e2,  e3)  frame  associated  with  orthotropy,  the  trans¬ 
formed  stress  tensor  of  Eq.  (10)  can  be  written  as 

Gmn  —  Lmnld^kl  =  ^mnkl^ klij  G ij  (20) 


where  K  is  the  fourth-order  deviatoric  projection  operator  whose 
components  with  respect  to  any  Cartesian  coordinate  system  are 

Km  =  l(8ikdji  +  8ndjk)-^pL-  (21) 

In  Voigt  notation,  K  can  be  written  as 
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Substituting  the  anisotropic  stress  potential  given  by  Eq.  (14)  into 
the  flow  rule  of  Eq.  (1)  and  applying  the  chain  rule  yields 


dtp  _  •  dtp  ddmn  da'rs  _  •  dcp  T  T/ 
dOij  ~Ad&mn  d(7'rs  dffij  ~  d&  LmnrsKr 


Hence,  in  the  anisotropic  case, 


(23) 


Brsijdi]  =  when  (p{a^d2,a 3;  rh,  a\)  =  0,  (24) 

where  the  fourth-order  tensor  B  =  L_1.  Indeed,  if  the  constraints  gi¬ 
ven  by  Eq.  (18)  are  enforced  then  L  is  invertible.  Note  that  in  the  iso¬ 
tropic  case, 


dy  =  4o when  <p(s, ,s2,s3;m, oT)  =  0,  (25) 

where  m  is  the  isotropic  version  of  m  given  by  Eq.  (16)  and  reduces 
to 


m  =  m|t_i  =  J  2  9/2 - .  (26) 

V  3k2  -2k +  3 

The  plastic  multiplier  rate  associated  with  the  isotropic  CPB06 
stress  potential,  Aiso,  is  expressed  in  terms  of  all  principal  values  of 
the  stress  deviator  (see  Cazacu  and  Stewart,  2009).  If  d/  ^  dn  ^  dIU 
are  the  principle  values  of  the  rate  of  deformation  tensor  d,  then 
the  expression  for  the  plastic  multiplier  rate  is  given  as 


■^iso  — 


(1  /  (3k2+10k+3\^  \  rl2  \  rl2 

WM)  y  t  3fc^-2k+3  J  dl  +  +  am 

1  Id 2  id2  I  /^3/<2-10k+3^2 

ip  Y  +  a"  +  \  3k2+2/f+3  ) 


if  d/  >  3k2-2k+3 

y/dijdij  ^  y/6(k2+3)(3k2+l)  ’ 

if  djn  ^  -(3k2+2k+3) 

y/<Wij  ^  y/ 6(/<2+3)(3k2+l)  ‘ 

(27) 


If  k  =  0  (no  tension-compression  asymmetry)  then  A  reduces  to  the 
von  Mises  effective  strain  rate  (A  =  v/(2/3)d  :  d).  Comparing  Eqs. 
(24)  and  (25),  it  follows  that  in  the  anisotropic  case  the  expression 
for  the  plastic  multiplier  rate  A  is  obtained  by  replacing  in  Eq.  (27) 
the  rate  of  deformation  tensor  d  with  b  =  B:d  and  the  constant  m 
with  m.  Thus,  the  anisotropic  plastic  multiplier  rate  is  given  as 
follows: 


jf  >  3k2-2k+3 

y/bijbij  ^  6(k2+3)(3k2+l)  ’ 

if  hi  ^  -(3k2+2k+3) 

y/ 6(/<2+3)(3k2+l), 

(28) 


where  bh  bn  and  bIU  are  the  ordered  principal  components  of  by  (for 
more  details  see  Cazacu  et  al.,  2010).  It  is  worth  noting  that  if  there 
is  no  tension-compression  asymmetry  in  the  matrix  material  (i.e., 
the  yield  strength  in  tension  is  equal  to  the  yield  in  compression) 
then  the  parameter  k  associated  with  strength  differential  effects 
is  automatically  zero.  Therefore,  the  anisotropic  strain  rate  potential 
given  by  Eq.  (28)  reduces  to  A  =  Vb  :  b/rh,  which  can  be  rewritten 
as  A  =  y/ (2/3)d  :  H  :  d  with  H  being  a  fourth-order  orthotropic  ten¬ 
sor  satisfying  Hiiki  =  0  with  i  =  1,  2,  3.  Thus,  if  k  =  0,  Eq.  (28)  reduces 
to  Hill’s  (1987)  orthotropic  strain  rate  potential. 

Note  that  the  scalar  b:b,  which  appears  in  Eq.  (28),  can  be  ex¬ 
pressed  as  follows: 


b:b  =  d:L:d 


(29) 


where  L  is  a  fourth-order  diagonal  tensor  in  the  coordinate  system 
associated  to  the  material  symmetry  (elt  e2,  e3)  and  is  expressible  in 
Voigt  notation  as 

~h  o  0  0  0  o' 


0  l2  0  0  0  0 
o  o  i3  o  o  o 
0  0  0  u  o  0 
o  o  o  o  i5  o 
0  0  0  0  0  }6 


(30) 


The  expressions  for  the  diagonal  components  of  L  in  terms  of  the 
components  of  B  =  L  1  are 


h  =  +  #i2  +  #i3  -  Ei  1E12  -  BnBn  -  E22E12  +  B22B23  -  E33E13 

+  B33B23  +  E12B13  —  E12B23  —  E13B23 

h  =  B\2  +  B222  +  B2 3  -  B22E12  -  B22B23  -  BuBu  +  BuBn  -  B33B2 3 

+  E33B13  +  B12E23  —  E23E13  —  B12B13 
h  =  B\ 3  +  B\3  +  B33  —  B33B13  —  B33B23  —  E11B13  +  B11B12  —  B22B23 
+  E22E12  +  E13B23  —  E12B13  —  B23E12 
U  =  B244 
Is  =  B55 
^6  =  B66. 

(31) 


4.  Development  of  the  anisotropic  macroscopic  yield  criterion 

An  approximate  yield  criterion  will  now  be  analytically  derived 
for  a  void-matrix  aggregate  containing  randomly-distributed 
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spherical  voids  and  with  an  incompressible  matrix  material  dis¬ 
playing  anisotropy  and  tension-compression  asymmetry.  The  ma¬ 
trix  yield  behavior  will  be  described  using  the  criterion  given  by 
Eq.  (14).  The  approach  that  will  be  used  in  the  derivation  is  the 
kinematic  homogenization  approach  presented  in  Section  2. 

Consider  as  an  RVE  a  hollow  sphere  with  inner  radius  a  and  out¬ 
er  radius  b  =  a/_1/3,  where  /is  the  void  volume  fraction  (also  called 
the  porosity).  The  void  is  considered  to  be  traction-free.  The  RVE  is 
subjected  to  axisymmetric  loading  such  that 

L  =  <g) ei  +  e2  0e2)  +  £33(e3  0e3),  (32) 

where  £  =  (a)a  denotes  the  macroscopic  stress  tensor  while  (ei,  e2, 
e3)  is  the  frame  associated  to  the  axes  of  orthotropy. 

To  obtain  the  upper-bound  of  the  overall  plastic  potential,  the 
classical  velocity  field  proposed  by  Rice  and  Tracey  (1969)  and 
Gurson  (1977)  will  be  used.  Thus,  the  local  velocity  field  v  in  the 
RVE  is  considered  to  be  of  the  form 

v  =  vu+vs,  (33) 


4A.  Overall  plastic  dissipation 

Let  W*(D)  denote  the  average  macroscopic  plastic  dissipation 
corresponding  to  the  velocity  field  given  by  Eq.  (35)  such  that 

W+(D)=1  [  w(d)dV  =  L  [  a'MV,  (39) 

v  Jq-co  *  Ja-(o 

where  d  is  given  by  Eqs.  (36)  and  (37),  V  =  ^nb3  is  the  volume  of  the 
RVE  being  considered,  co  represents  the  void  volume  and  A  is  the 
plastic  multiplier  rate  associated  to  the  anisotropic  CPB06  yield  cri¬ 
terion  (see  Eq.  (28)).  The  approach  of  Leblond  (2003)  will  be  fol¬ 
lowed  to  obtain  a  new  upper-bound  estimate  of  the  macroscopic 
yield  locus  by  applying  the  Cauchy-Schwartz  inequality  as  follows: 

W+(D)=1^47ir2(w(d))S(r)dr 

<  W++(D)  =  i  J\nr2  [<w2(d))S(r)]1/2dr,  (40) 


where  vv  is  associated  with  volumetric  expansion  and  Vs  is  associ¬ 
ated  with  shape  changes.  Imposing  uniform  rate  of  deformation 
boundary  conditions  and  matrix  incompressibility  yields 

v(x  =  ber)  =  D  x  and  divv  =  0,  (34) 

where  x  is  the  Cartesian  position  vector  that  denotes  the  current  po¬ 
sition  in  the  RVE  and  er  is  the  radial  unit  vector. 

From  Eqs.  (33)  and  (34),  it  follows  that 


Vs  =  D  x, 


where  r  is  the  radial  coordinate,  Dm  =  ±Dkk  is  the  mean  part  of  D  and 
D'  =  D  -  Dm I  is  the  deviatoric  part  of  D  with  I  being  the  second-or¬ 
der  identity  tensor.  In  Eq.  (35),  both  Cartesian  and  spherical  coordi¬ 
nates  have  been  used.  The  local  rate  of  deformation  tensor 
d  =  \  ( Vv  +  VTv)  corresponding  to  the  velocity  field  v  given  by 
Eq.  (35)  is 

d  =  dv  +  ds,  (36) 

where 


ds  =  D , 


with  d  =  (-2er  0  er  +  e0  0  e0  +  e^,  0  e^).  Hence,  the  principal  values 
of  the  local  rate  of  deformation  tensor  d  are  as  follows: 


where  S(r)  is  the  spherical  surface  of  radius  r.  In  Eq.  (40),  the  follow¬ 
ing  notation  was  used  for  averaging  over  S(r): 

-i  phi  pH 

(*>s<r)=4 nj0  l  xsin9ded(l>-  (41) 

When  the  matrix  behavior  is  described  by  the  von  Mises  yield 
criterion  (L  =  l,  k  =  0),  Wf+(D)  can  be  evaluated  analytically  (see 
Gurson,  1977;  Leblond,  2003).  In  the  current  derivation,  fresh  dif¬ 
ficulties  are  encountered  when  estimating  the  overall  local  plastic 
dissipation  w(d).  These  difficulties  are  due  to  the  tension-com¬ 
pression  asymmetry  and  anisotropy  of  the  matrix  material  and 
are  associated  with  the  fact  that  the  CPB06  plastic  multiplier  rate 
has  multiple  branches  (see  Eq.  (28)).  The  expression  for  the 
CPB06  plastic  multiplier  rate  depends  on  the  sign  and  ordering  of 
the  principal  values  of  the  transformed  rate  of  deformation  tensor 
b  =  B:d. 

For  the  case  of  an  isotropic  matrix  exhibiting  tension-compres¬ 
sion  asymmetry  (L  =  1  but  k  ^  0),  the  expressions  for  the  local  plas¬ 
tic  dissipation  are  provided  in  Cazacu  and  Stewart  (2009).  In  order 
to  obtain  an  analytic  expression  of  the  integral  representing  the 
overall  plastic  dissipation,  it  was  assumed  that  in  the  expressions 
of  the  local  plastic  dissipation  (see,  for  example,  Eqs.  (24)-(30)  of 
Cazacu  and  Stewart,  2009),  the  cross  term  DmD'^^  could  be 
neglected.  Note  that  the  upper-bound  character  of  the  criterion  is 
retained  for  both  the  purely  deviatoric  or  purely  hydrostatic  cases 
(since  the  cross  term  is  zero  for  both  loading  conditions).  The  valid¬ 
ity  of  the  approximation  for  general  loadings  was  assessed  in 
Cazacu  and  Stewart  (2009)  by  conducting  finite-element  cell 
calculations.  Using  this  approximation,  the  isotropic  local  plastic 
dissipation  becomes  (see  Eqs.  (19),  (27)  and  (38)): 

w(d)  =  sjndfdij,  (42) 

where 


d2  =  - 


1 

2 


D'n 


+ 


2D’UD„ 


cos  2  0+  (D),)2, 


d3 


-0 


+  2DnDm 


cos  2 9  3-  (D'n)2. 


(38) 


if  <  0, 

if  Zm  >  0, 


(43) 


and  Zm  is  the  applied  mean  stress.  Note  that  for  the  cases  of  purely 
hydrostatic  and  purely  deviatoric  loading,  the  approximation  given 
by  Eq.  (42)  coincides  with  the  exact  value  of  the  local  plastic  dissi¬ 
pation  for  the  velocity  field  given  by  Eq.  (35). 

In  the  anisotropic  case,  the  local  plastic  dissipation  is  obtained 
by  replacing  in  Eq.  (42)  the  local  rate  of  deformation  tensor  d  with 
the  transformed  rate  of  deformation  tensor  b  =  B:d  and  the  con¬ 
stant  n  with  its  anisotropic  equivalent  h.  Thus,  the  following 
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approximate  form  of  the  plastic  multiplier  rate  will  be  used  when 
calculating  the  overall  plastic  dissipation  W++  (see  Eq.  (40)): 

A2  =  hbijbij 


such  that 

w(d)  =  g\ 

\J  nbjjbjj, 

where 

(  —  1 

<  3  A 

Hx, 

\3k2-2k+3  J 
f  3  A 

[  m2  \ 

\3k2  +2k+3  J 

(46) 


4.2.  Macroscopic  stresses 

Let  Ie  denote  the  macroscopic  effective  stress  associated  to  the 
CPB06  anisotropic  potential  as 


IP  =  m 


E(i^i-^,)2 


(47) 


where  E  =  L  :  E'  is  the  transformed  macroscopic  stress  tensor.  The 
effective  strain  rate  associated  with  the  CPB06  stress  potential  will 
be  denoted  by  De.  A  rigorous  proof  that  De  is  a  work-conjugate  mea¬ 
sure  of  Ie  (i.e.,  that  IeDe  =  E' :  D')  was  provided  in  Cazacu  et  al. 
(2010).  According  to  Eqs.  (29)  and  (44) 


De  =  A(D')  =  Jn(D':L:D'). 


(48) 


Let  Ie  =  ^(3/2 )I'ijI'ij  and  De  =  ^(2/3 )D'ijD'ij  denote  the  von 
Mises  effective  stress  and  the  von  Mises  effective  strain  rate, 
respectively.  For  an  orthotropic  material  under  axisymmetric  load¬ 
ing  conditions  (see  Eq.  (35)),  the  von  Mises  effective  stress  is  re¬ 
lated  to  the  CPB06  effective  stress  as 


^ e  —  S^e 


(49) 


where  g  is  a  constant  (see  Eq.  (75))  expressible  in  terms  of  the 
strength  differential  parameter  k  and  the  components  of  the  aniso¬ 
tropic  tensor  L.  Since  De  is  a  work-conjugate  measure  of  Ie  and  Ie  is 
a  work-conjugate  measure  of  De,  it  follows  that 


De  =  -De. 
g 


(50) 


The  macroscopic  stress  tensor  associated  with  the  upper-bound 
estimate  of  the  plastic  dissipation  VlA+(D)  is  given  as 


dW++ 

x=—iDm,De). 

Employing  the  chain  rule  yields 

dW++  dDm  dW++  dDe 


L  =  '  dDm  dD 

such  that 

_  1  dW++ 

2*  m  —  ' 


dDe  dD 


3  dDm 


and 


v/  2  D  dW++ 
“3  De  dDe  ■ 
It  follows  that 

dW++ 


Ie  = 


dDe 


(51) 


(52) 


(53) 


(54) 


(55) 


and 

dW++  1  dW++ 


(44) 

dDe 

g  dDe 

such  that 

(45) 

y  — 

dW++ 

Zje  — 

dDe 

(56) 


(57) 


5.  Anisotropic  criterion 

According  to  Eq.  (37),  the  local  rate  of  deformation  can  be 
decomposed  into  a  volumetric  part  dv  and  a  deviatoric  part 
ds  =  D'  =  constant.  Thus,  the  approximate  expression  for  the  macro¬ 
scopic  plastic  dissipation  can  be  written  as 


dr 

1  1/2 


W++ =  $- £  4*r2[(X2)m] 

=  y  l  4ra-2[n(b:b)S(r)j 
=  if  [  47tr2[ft<d:L:d)s(r)]''"* 

=  y£ [n^  :  L  :  dv  +  ds  :  L  :  ds  +  2dv  :  L  :  ds) 


dr 
"1 1/2 


1  1/2 


S(r)  J 


Since 

(d:LdEn=D:L:D 
and  noting  that 

(dl,:L:dS)s(r)  =  (dV)s(r):L:D' 
with 

(d)S(r)  =  ( [-2er  ®  er  +  e„  ®  e9  +  e*  ® 


)s(n 


dr. 

(58) 

(59) 

(60) 

(61) 


the  following  expression  is  obtained  for  the  estimated  macroscopic 
plastic  dissipation: 


W++  =- 


-  f  4 [h^V  ■■ 1  ■■  dV>s(r)  +  nD' :  L  :  O'] 1/2 dr 


r2\lDl^h2  +  DMr 


M  (b 

b3  Ja 

^|Dm|h^V/yC 


=  G 

such  that 

W++=a\\Dm\h 
where  u  =  b3jr3, 

y  =  mh 

De 

and 

h2  =  n(d  :  L  :  d>S(r)  =  n 


+ 


1  du 

y2  u2 


i  yj\  +  y2u2 

sinh  (yu)  -  — - - — 

w  7  yu 


i  // 


f+i2  +  i3)+^(i4  +  h  +  ie) 


(62) 


(63) 


(64) 


(65) 


(see  Appendix  A  for  details  regarding  the  calculation  of  h).  Note  that 
the  parameter  h  depends  on  the  anisotropy  coefficients  (see  Eq. 
(15))  as  well  as  the  parameter  h  given  by  Eq.  (46). 
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Using  Eqs.  (53)  and  (57)  along  with  Eq.  (63)  yields 

’++  U  i  \  if 

2  cr[sgn(Dm)  [sinh-1  (yu)J  ^ 


1  dW 

^ m  — 


3  dDm  3 
=  |c[sgn(Dm) 


sinh 


1  (j)  -  sinh  1  (y) 


and 


dW+ 


dDe 


=  <7 


x/l  +  y2u2 


i  // 


= 


V1  +p 


(66) 


(67) 


Finally,  eliminating  y  from  the  previous  two  equations  yields  the 
following  expression  for  the  macroscopic  yield  criterion  incorporat¬ 
ing  anisotropy: 


U 


0  = 


mJEt,  l^.-l  -kZ, 


Y 


(68) 

where  t  =  L  :  is  the  uniaxial  tensile  yield  strength  in  the  1- 

direction  of  orthotropy,/is  the  void  volume  fraction,  Zm  is  the  mac¬ 
roscopic  mean  stress,  h  is  the  hydrostatic  factor  given  by  Eq.  (65) 
and  the  parameter  m  is  given  by  Eq.  (16).  These  key  equations  are 
repeated  here  for  the  readers  convenience: 


/=(i-/)D^  =  (i-/)ig- 

-‘'-•KrWr)-  <69> 

Eq.  (68)  represents  the  approximate  yield  criterion  for  a  porous 
material  with  the  matrix  incompressible,  orthotropic  and  display¬ 
ing  tension-compression  asymmetry.  Although  the  expression  of 
the  developed  criterion  is  similar  to  that  of  Benzerga  and  Besson 
(2001),  there  are  distinct  differences.  First,  Hill’s  (1948)  effective 
stress  is  replaced  with  the  CPB06  effective  stress  Ze,  which  involves 
the  stress  deviator  27 ,  the  anisotropy  coefficients  Ly  and  the 
strength-differential  parameter  k.  Secondly,  the  hydrostatic 
parameter  h  depends  on  the  anisotropy  coefficients  Ly,  on  the 
strength-differential  parameter  k  and  on  the  sign  of  the  applied 
mean  stress  Zm. 

In  contrast  with  the  criteria  of  Gurson  (1977)  and  that  of  Benze¬ 
rga  and  Besson  (2001 ),  the  yield  locus  given  by  Eq.  (68)  is  no  longer 
symmetric  with  respect  to  the  axis  Zm  =  0.  According  to  the  devel¬ 
oped  criterion,  for  tensile  hydrostatic  loading,  yielding  occurs  in 
the  porous  aggregate  when  Em  =  pj  with 


Py  = 


7h  +2/3  +  61, 


ln(f), 


(70) 


Y  15m2  \3k2  +  2k  +  3 
while  for  compressive  loading  yielding  occurs  when  =  Py  with 


/  1  (16) 

2  1 

( 7/i  +  2/3  +  6/4\ 

\l(\0,\-k0,)2  +  (\02\-k02)2  +  (m-k03)2  Py=^a 

15m2  1 

(  3k2  -  2k  +  3  J 

m  = 


with 

.  2  1  1 

V]  =  jMi  -  -L12  -  JM3, 

,  2  1  1 

(l}2  =  ^M2  —  ^L22  ~  ^  ^23  • 

,  2  1  1 

^3  =  13  —  ^  f<23  —  ^  ^33 


I  ln(f). 


(71) 


(17) 


and 

h2  =  n(d  :  L  :  d)S(r)  =  n 

with 


i(i,+i,+ 


--(I4  +  I5  + 


k) 


(65) 


If  there  is  no  tension-compression  asymmetry  in  the  matrix  (k  =  0), 
pf  =  \py  |.  The  limiting  pressures  p+  and  p~  correspond  to  the  exact 
solution  of  a  hollow  sphere  loaded  hydrostatically  only  if  the  matrix 
is  isotropic  (L  =  1)  since,  in  this  case,  the  trial  velocity  field  given  by 
Eq.  (35)  is  the  only  admissible  velocity  field. 

If  the  matrix  is  orthotropic  and  has  no  tension-compression 
asymmetry,  Eq.  (68)  reduces  to  that  of  Benzerga  and  Besson 
(2001).  Indeed,  for  k  =  0,  h  =  2/3  and  L  =  H  (i.e.,  the  matrix 
behavior  is  described  by  Hill’s  1948  yield  criterion),  the  hydrostatic 
factor  h  given  by  Eq.  (65)  reduces  to 


h  =  \l^(ii+h  +  h)+$(u+i5  +  k), 


15 


(72) 


h  =  +  £>12  +  fi?3  -  BnBn  -  #11813  -  #22^12  +  B22B23  -  £33^13 

+  B33B23  +  BnBu  —  BuB23  —  #13^23 
h  =  B2U  +  #22  +  B\3  —  B22B\2  —  B22B23  —  ^11^12  +  B-[  1 B13  —  B33B23 

+  ^33^13  +  BuB23  —  #23^13  ~  ^12^13 

l3  =  B\3  +  833  +  ^33  —  B33B-13  —  B33B23  —  #11  #13 
+  BnB\2  +  B\3B23  —  BUBu  —  B23B\2 
l  =  Bl 


which  is  the  hydrostatic  factor  corresponding  to  a  Hill  (1948)  ma¬ 
trix  (see  Benzerga  and  Besson,  2001).  Note  that  when  the  matrix 
is  isotropic  and  displays  tension-compression  asymmetry  (i.e., 
L  =  I  and  k  Y  0), 


B11B12  —  B22B23 


h  = 


2  if  I’m  <  0, 

2ym  if^>° 


(73) 


h  =  Bj55 
fe  =  B66, 

(31) 

where  B  =  L_1. 

Note  that  the  developed  criterion  (Eq.  (68))  accounts  for  the 
combined  effects  of  matrix  anisotropy  and  tension-compression 
asymmetry  on  both  the  yielding  of  the  porous  aggregate  and  the 
porosity  evolution.  These  effects  are  captured  through  the 
strength-differential  parameter  k  and  the  anisotropy  coefficients 
(the  components  of  the  fourth-order  tensor  L)  involved  in  the 
expressions  of  the  transformed  stress  tensor  d  as  well  as  m  and 
h.  Indeed,  the  evolution  of  the  porosity  is  given  as 


such  that  the  developed  criterion  (Eq.  (68))  reduces  to  that  pro¬ 
posed  in  Cazacu  and  Stewart  (2009).  Furthermore,  for  a  von  Mises 
matrix  (L  =  I  and  k  =  0)  the  developed  criterion  reduces  to  Gurson’s 
(1977)  yield  criterion  for  spherical  voids.  In  the  next  section,  the 
capabilities  of  the  developed  model  will  be  further  illustrated  and 
comparisons  with  finite  element  calculations  will  be  made  to  vali¬ 
date  the  predicted  yield  behavior. 

6.  Validation  of  the  developed  anisotropic  criterion 

This  section  focuses  on  the  assessment  of  the  developed  aniso¬ 
tropic  criterion  given  by  Eq.  (68).  The  model  is  analyzed  in  detail 
for  some  specific  materials  in  Section  6.1,  axisymmetric  unit  cell  fi¬ 
nite  element  (FE)  calculations  are  presented  in  Section  6.2  and 
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comparisons  between  the  finite  element  results  and  the  materials 
analyzed  previously  are  given  in  Section  6.3. 

6.1.  Model  assessment 

The  following  analysis  will  focus  on  the  case  of  transversely  iso¬ 
tropic  materials.  The  choice  of  a  simplified  material  orthotropy 
representation  (i.e.,  the  plane  of  symmetry  eie2  is  isotropic)  is  con¬ 
sidered  because  the  discussion  for  general  orthotropy  would  in¬ 
volve  a  very  high  number  of  materials  corresponding  to  different 
values  for  the  anisotropy  coefficients  Ly  (i.e.  various  orderings  of 
the  tensile  and  compressive  yield  stresses  corresponding  to  the 
directions  ei,  e2,  e3).  Transverse  isotropy  is  common  in  clock-rolled 
hexagonal  metal  sheets,  where  multiple  passes  ensure  in-plane 
isotropy  (see,  for  example,  the  data  on  high-purity  zirconium  re¬ 
ported  in  Tome  et  al.,  2001 ;  Nixon,  2008).  For  a  given  loading  path, 
the  parameters  that  are  varied  are  the  anisotropy  coefficients  Ly 
and  the  tension-compression  parameter  k.  Thus,  the  yield  surfaces 
of  nine  materials  containing  randomly  oriented  spherical  voids  will 
be  examined: 

1.  Isotropic  materials  for  which  the  matrix  is  characterized  by 
L  =  1  and  either  displays  tension-compression  asymmetry 
(such  that  k  ^  0)  or  does  not  (/<  =  0).  These  materials  are  called 
materials  type  A. 

2.  Transversely  isotropic  materials  for  which  the  matrix  has  a 

weaker  in-plane  yield  strength  than  through-thickness  yield 
strength  (see  Fig.  2)  with  either  tension-compression  asymme¬ 
try  (/<  0)  or  not  (/<  =  0).  These  materials  are  called  materials 

type  B. 


3.  Transversely  isotropic  materials  with  matrix  displaying  larger 
in-plane  yield  strength  than  through-thickness  yield  strength 
(see  Fig.  2)  with  either  tension-compression  asymmetry 
(/<  0)  or  not  (k  =  0).  These  materials  are  called  materials 

type  C. 

The  isotropic  case  (material  type  A)  is  taken  as  a  reference.  Note 
that  calculations  for  material  A  with  k  =  0  corresponds  to  a  von 
Mises  matrix;  hence,  the  FE  results  obtained  in  this  case  will  be 
compared  to  the  Gurson-Tvergaard-Needleman  (GTN)  yield  locus. 
Likewise,  material  type  A  with  k  ^  0  corresponds  to  the  isotropic 
CPB06  matrix  and  FE  results  will  be  compared  to  the  Cazacu  and 
Stewart  (2009)  yield  locus.  In  Table  1  are  given  the  numerical  val¬ 
ues  of  the  anisotropy  coefficients  for  each  material.  These  values 
were  chosen  such  that  all  materials  have  the  same  tensile  yield 
strength  in  the  plane  of  symmetry.  As  an  example,  in  Fig.  2  are  rep¬ 
resented  the  projection  of  the  yield  loci  given  by  Eq.  (68)  corre¬ 
sponding  to/=  0  (the  void-free  state)  for  all  these  materials. 

As  already  mentioned,  if  the  anisotropic  matrix  does  not  display 
tension-compression  asymmetry  (/<  =  0),  then  the  developed  crite¬ 
rion  of  Eq.  (68)  reduces  to  that  proposed  by  Benzerga  and  Besson 
(2001)  for  spherical  void  geometry.  In  Benzerga  and  Besson 
(2001),  through  comparison  with  axisymmetric  cell  calculations, 
the  authors  have  assessed  the  accuracy  of  the  description  of  the  ef¬ 
fects  of  plastic  anisotropy  on  void  growth.  However,  no  discussion 
of  the  accuracy  of  the  yielding  description  has  been  reported. 

In  this  paper,  the  numerical  values  for  the  anisotropy  coeffi¬ 
cients  of  the  transversely  isotropic  materials  B  and  C  with  no  ten¬ 
sion-compression  asymmetry  were  taken  to  coincide  with  the  set 
of  anisotropy  parameters  considered  by  Benzerga  and  Besson 
(2001),  which  correspond  to  a  zircaloy  sheet  and  thin  aluminum 


Fig.  2.  Plane  stress  yield  loci  for  void-free  materials  A  (isotropic),  B  and  C  (transversely  isotropic)  according  to  the  CPB06  yield  criterion,  x  is  an  in-plane  direction  with  z  being 
the  through-thickness  direction,  (a)  No  strength  differential  ( k  =  0).  (b)  Tension-compression  asymmetry  with  the  yield  strengths  in  tension  less  than  the  yield  strengths  in 
compression  (/<  =  -0.3098).  (c)  Tension-compression  asymmetry  with  the  yield  strengths  in  tension  greater  than  the  yield  strengths  in  compression  ( k  =  0.3098). 
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Table  1 

Numerical  values  of  the  anisotropy  coefficients  for  the  materials  A,  B,  and  C. 


CPB06  parameters 

Material  A 

Material  B 

Material  C 

f-n  =  L22 

1.000 

1.054 

0.963 

L33 

1.000 

0.850 

1.064 

f-13  =  f-23  =  2O  -  T33) 

0.000 

0.075 

-0.032 

Lu=Ul+L33)-Lu 

0.000 

-0.129 

0.069 

L44 

1.000 

0.775 

1.817 

sheet,  respectively.  This  choice  of  anisotropic  coefficients  serves 
two  purposes:  (a)  it  allows  partial  verification  of  the  FE  implemen¬ 
tation  of  the  model  (Eq.  (68))  and  (b)  it  provides  FE  data  for  assess¬ 
ing  the  analytic  yield  loci  according  to  Benzerga  and  Besson  (2001) 
for  k  =  0  and  according  to  Eq.  (68)  for  k  ^  0. 

For  axisymmetric  loadings,  the  effective  stress  according  to  the 
CPB06  anisotropic  criterion  is  related  to  the  von  Mises  effective 
stress  by  the  relation 


where 


g  = 


_  i 

m|(P3|y/2(3k2+2k+3) 

.  f 

m\03\^/2(3k2-2k+3) 


for  X\\  >  2^33, 
for  Zu  <  ^33- 


(75) 


In  Fig.  3  are  presented  the  projections  of  the  analytical  yield  loci  gi¬ 
ven  by  Eq.  (68)  in  the  plane  (Ze,  Xm)  for  a  transversely  isotropic 
material  of  type  B  with  either  yield  strengths  in  tension  greater  than 
in  compression  (/<  =  0.3098)  or  yield  strengths  in  tension  less  than  in 
compression  (/<=  -0.3098). 

Note  that  the  effect  of  the  sign  of  the  third  invariant  of  the  stress 
deviator  is  well  captured.  Specifically,  for  purely  deviatoric  loading, 
if  X33  >  Zn  then  the  material  yields  at 


Ze=~ 


<4(1  -/) 


m|$3|y,2(3k2  -2k  +  3)’ 


(76) 


=  |^11  “  ^33  |  =  gie. 


(74) 


whereas  for  Z33  < 


(a) 


w" 

I 


E-h  _ 

b 


2Zn  +S33 


2^ll  +^33 

3  erf 


Fig.  3.  Anisotropic  yield  surfaces  for  a  matrix  material  containing  a  spherical  void  of  porosity/=  0.01  and  with  the  through-thickness  yield  strength  greater  than  the  in-plane 
yield  strength  (i.e.,  material  B ).  (a)  Tension-compression  asymmetry  with  the  yield  strengths  in  tension  less  than  the  yield  strengths  in  compression  (/<  =  -0.3098).  (b) 
Tension-compression  asymmetry  with  the  yield  strengths  in  tension  greater  than  the  yield  strengths  in  compression  (k  =  0.3098). 
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v  m  -/) 

m|4>3|v/2(3fc2+2fc  +  3)’ 
where,  for  transverse  isotropy, 

^>3  =  g(l  —  3  L33). 


(77) 


(78) 


Note  also  that  the  yield  locus  is  no  longer  symmetric  with  the  ver¬ 
tical  axis  Zm  =  0  since  h  depends  on  the  sign  of  the  mean  stress  (Eq. 
(65)).  For  transverse  isotropy,  h  can  be  written  as 


h  =  \/-n(7U 


-  2U  +  6L 


(79) 


The  effect  of  the  magnitude  of  the  third  invariant  of  the  stress 
deviator  can  be  best  assessed  through  examination  of  the  projec¬ 
tions  of  the  yield  locus  (Eq.  (68))  in  the  7i-plane  (the  plane  normal 
to  the  hydrostatic  axis:  (ei  +  e2  +  e3)/\/3).  Contours  corresponding 
to  a  porosity  of  /=  0.01  and  fixed  values  of  the  mean  stress 
Zm  =  0, 0.75py,  and  0.9 pj  are  plotted  in  Figs.  4-6  for  the  six  mate¬ 
rials  considered. 

Fig.  4(a)-(c)  shows  the  7i-plane  representation  of  the  yield  loci 
for  a  porous  aggregate  with  an  isotropic  matrix  (material  type  A) 


displaying  tension-compression  asymmetry  with  k  =  0.3098  and 
k  =  -0.3098,  along  with  the  von  Mises  matrix  {k  =  0)  for  compari¬ 
son.  The  effect  of  the  third  invariant  of  the  stress  deviator  is  evi¬ 
dent  in  the  change  in  shape  of  the  yield  loci  from  circles  to 
rounded  triangles.  The  combined  effects  of  anisotropy  and  ten¬ 
sion-compression  asymmetry  are  evident  from  Fig.  5  (material 
type  B)  and  Fig.  6  (material  type  C).  A  very  drastic  departure  of 
the  yield  loci  from  the  ellipse  (/<  =  0)  is  observed  for  both  material 
types  B  and  C. 

Fig.  7  shows  the  7i-plane  representation  of  the  ductile  criterion 
given  by  Eq.  (68).  Fig.  7(a)  illustrates  the  7r-plane  representation  for 
a  von  Mises  material  (material  type  A  with  k  =  0  and  /=  0).  Figs. 
7(b)-(d)  shows  representations  with  non-zero  porosity  (/=  0.01 ) 
for  materials  of  type  B  with  no  strength-differential  ( k  =  0),  yield 
strength  in  tension  greater  than  in  compression  (k  >  0)  and  yield 
strength  in  tension  less  than  in  compression  (/<  <  0),  respectively. 
Notice  that  if  there  are  no  voids  present  (see  Fig.  7(a)),  the  response 
is  independent  of  the  hydrostatic  pressure;  this  is  because  the  ma¬ 
trix  material  is  incompressible.  When  voids  are  present,  as  in 
Fig.  7(b)-(d),  the  yield  surface  shrinks  in  the  deviatoric  plane  as 
the  pressure  increases  toward  either  the  tensile  or  compressive 
hydrostatic  limit  pressure  (i.e.,  toward  either  vertex).  Fig.  7  implies 
that  the  anisotropy  and  tension-compression  asymmetry  of  the 
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Fig.  4.  Representation  in  the  deviatoric  plane  of  the  ductile  yield  criterion  given  by  Eq.  (68)  for  an  isotropic  material  (material  A  in  Table  1 ).  p;  is  the  tensile  hydrostatic  yield 
pressure,  (a)  No  strength  differential  (k  =  0).  (b)  Tension-compression  asymmetry  with  the  yield  strengths  in  tension  less  than  the  yield  strengths  in  compression 
( k  =  -0.3098).  (c)  Tension-compression  asymmetry  with  the  yield  strengths  in  tension  greater  than  the  yield  strengths  in  compression  (/<  =  0.3098). 
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Fig.  5.  Representation  in  the  deviatoric  plane  of  the  ductile  yield  criterion  given  by  Eq.  (68)  for  a  material  with  the  through-thickness  yield  strength  greater  than  the  in-plane 
yield  strength  (material  B  in  Table  1).  pj  is  the  tensile  hydrostatic  yield  pressure,  (a)  No  strength  differential  (/<  =  0).  (b)  Tension-compression  asymmetry  with  the  yield 
strengths  in  tension  less  than  the  yield  strengths  in  compression  ( k  =  -0.3098).  (c)  Tension-compression  asymmetry  with  the  yield  strengths  in  tension  greater  than  the  yield 
strengths  in  compression  (/<  =  0.3098). 


matrix  influences  the  shape  of  the  yield  locus  in  the  deviatoric 
plane  while  the  presence  of  voids  in  the  aggregate  leads  to  a 
decrease  in  the  size  of  the  yield  locus  with  increasing  pressure 
(positive  or  negative).  Fig.  7(b)-(d)  show  yield  surfaces  corre¬ 
sponding  to  a  fixed  value  of  the  porosity,  /=0.01;  varying  the 
porosity  will  simply  change  the  size  of  the  yield  locus  (the  yield  lo¬ 
cus  decreases  in  size  with  increasing  porosity  for  a  fixed  pressure). 
The  shape  of  the  yield  locus  is  governed  by  the  matrix  anisotropy 
and  tension-compression  asymmetry. 


(1997)  and  Cazacu  and  Stewart  (2009)  in  the  sense  that  the  contin¬ 
uum  is  considered  to  consist  of  a  periodic  assemblage  of  hexagonal 
cylindrical  unit  cells  which  are  approximated  by  right  circular  cyl¬ 
inders.  Due  to  symmetry,  only  one  quarter  of  the  unit  cell  for  the 
anisotropic  material  is  shown  in  Fig.  8.  In  the  figure,  theX3-axis  de¬ 
notes  the  through-thickness  direction  while  the  Xi-axis  is  in  the 
plane  of  isotropy.  Every  cell  of  initial  length  L0  and  radius  R0  con¬ 
tains  a  spherical  void  of  initial  radius  a0  and  is  subject  to  homoge¬ 
neous  radial  and  axial  displacements.  In  other  words,  the  boundary 
conditions  for  this  unit  cell  are  as  follows: 


6.2.  Validation  through  finite  element  cell  calculations 

The  assumption  of  in-plane  isotropy  allows  axisymmetric  calcu¬ 
lations  to  be  performed  and  the  results  can  be  used  to  examine 
directional  effects  and  the  influence  of  tension-compression  asym¬ 
metry  on  yielding.  Thus,  loading  paths  of  the  type  (Zn,  In,  Z33) 
will  be  considered,  with  e3  being  normal  to  the  plane  of  isotropy. 
Tension-compression  asymmetry  in  the  matrix  due  to  the  com¬ 
bined  effects  of  the  matrix  asymmetry  and  the  presence  of  the 
voids  will  be  assessed  by  carrying  out  calculations  for  both  com¬ 
pressive  (Zm  <  0)  and  tensile  (Zm  >  0)  axisymmetric  loadings. 

The  axisymmetric  unit  cell  calculations  that  were  carried  out 
are  similar  to  those  of  Koplik  and  Needleman  (1988),  Ristinmaa 


Ui  =  ur2  =  ur3  =  0  for  =  0, 
u3  =  wc\  =  ur2  =  0  for  X3  =  0, 
ut  =  Ui  for  Xi  =  R , 
u3  =  U3  for  X3=L/ 2, 


(80) 


where  tq  are  the  displacements  in  the  X/  directions  and  urz  are  the 
rotations  about  the  X,  axes.  The  void  is  considered  to  be  traction 
free.  Thus,  for  this  axisymmetric  unit  cell,  the  initial  porosity  f0  is 
defined  as 


fo  = 


4a30 


(81) 
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Fig.  6.  Representation  in  the  deviatoric  plane  of  the  ductile  yield  criterion  given  by  Eq.  (68)  for  a  material  with  the  through-thickness  yield  strength  less  than  the  in-plane 
yield  strength  (material  C  in  Table  1).  pj  is  the  tensile  hydrostatic  yield  pressure,  (a)  No  strength  differential  (/<  =  0).  (b)  Tension-compression  asymmetry  with  the  yield 
strengths  in  tension  less  than  the  yield  strengths  in  compression  (k  =  -0.3098).  (c)  Tension-compression  asymmetry  with  the  yield  strengths  in  tension  greater  than  the  yield 
strengths  in  compression  (k  =  0.3098). 


The  macroscopic  principal  strains,  the  effective  von  Mises  mac¬ 
roscopic  strain  Ee  and  the  strain  triaxiality  TE  are  defined  as 
follows: 


Eii  =  In 


E33  =  In 


=  In 


=  In 


^o  +  2LT3j 


(82) 


^  =  a/^e;. 


(83) 


where  F3  is  the  total  radial  force  at  =  R  and  F3  is  the  total  axial 
force  at  X3  =  1/2.  Due  to  axisymmetric  loading,  J3  =  -(2/27)(Iu- 
X33)3;  thus,  only  the  effect  of  the  sign  of  the  third  invariant  can 
be  examined  in  the  current  analysis. 

The  analytical  yield  criterion  of  Eq.  (68)  can  be  modified  to  in¬ 
clude  additional  parameters,  as  was  done  by  Tvergaard  (1981) 
and  Tvergaard  and  Needleman  (1984)  in  the  case  of  Gurson’s 
(1977)  yield  criterion  as  follows: 

+  2<Ji/cosh  -  1  -  <hf2,  (86) 


and 


t  _  ^kk  _  2En  +E33 

£_3i;_2|£11  —  £33  r 


(84) 


where  Ekk  is  the  trace  of  the  macroscopic  strain  tensor  and  E'y  are  the 
components  of  the  macroscopic  deviatoric  strain  tensor. 

The  macroscopic  stress  £  is  also  axisymmetric  such  that 


27n 


(85) 


where  values  of  q3  =  4/e,  q2  =  1  and  q3  =  q\  will  be  used  in  the  fol¬ 
lowing  (these  values  were  suggested  in  Leblond  and  Perrin,  1990, 
based  on  an  analysis  of  the  effects  of  a  second  population  of  small 
voids  on  the  growth  of  a  large  void).  In  the  next  section,  the  analyt¬ 
ical  yield  surfaces  given  by  Eq.  (86)  will  be  compared  to  data  from 
finite  element  unit  cell  calculations. 

6.3.  Results 

All  calculations  were  performed  assuming  elastic,  ideal-plastic 
behavior  for  the  matrix  with  the  plastic  potential  given  by  the 
CPB06  anisotropic  criterion  (see  Eq.  (47)).  A  user  material 
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Fig.  7.  Representation  in  the  deviatoric  plane  of  the  ductile  yield  criterion  given  by 
Eq.  (68)  for  a  material  with  the  through-thickness  yield  strength  greater  than  the 
in-plane  yield  strength  (material  type  B  in  Table  1).  (a)  No  tension-compression 
asymmetry  (/<  =  0)  and  no  voids  (/=  0).  (b)  No  tension-compression  asymmetry 
with/=  0.01.  (c)  Yield  strengths  in  tension  greater  than  in  compression  ( k  =  0.3098) 
with/=  0.01.  (d)  Yield  strengths  in  tension  less  than  in  compression  (/<=  -0.3098) 
with/=  0.01. 


x3 


L/2 


*  X, 


Fig.  8.  Quarter  section  of  the  unit  cell  for  transversely  isotropic  materials 
containing  spherical  voids;  axis  X3  is  the  axis  of  rotational  symmetry  with  the 
Xi-X2  plane  being  the  plane  of  symmetry. 


subroutine  was  written  for  the  commercial  FE  program  ABAQUS 
(see  Abaqus,  2008)  using  a  fully  implicit  return  mapping  algorithm. 
A  mesh  refinement  study  was  performed  using  four-node  quadri¬ 
lateral  elements  (with  reduced  integration  and  ABAQUS’  enhanced 
hourglass  control)  to  ensure  solution  convergence.  Fig.  9  shows  the 
mesh  used  for  all  tests.  The  matrix  elastic  properties  are  E/o\  = 
800  and  v  =  0.32,  where  E  is  the  Young’s  modulus  and  v  is  the  Pois¬ 
son’s  ratio.  For  ease  of  comparison,  the  tensile  yield  strength  in  the 


Fig.  9.  /o  =  0.01  axisymmetric  finite  element  mesh  for  the  unit  cell. 


ei  direction,  o\,  was  considered  to  be  the  same  for  all  material 
types  A,  B,  and  C. 

The  model  implementation  was  first  verified  by  performing  cal¬ 
culations  corresponding  to  k  =  0  for  material  type  A  (isotropic)  and 
transversely  isotropic  (B  and  C),  which  correspond  to  a  von  Mises 
matrix  (for  material  type  A)  and  Hill  matrices  (for  material  types 
B  and  C),  and  comparing  the  results  to  those  obtained  using  the 
ABAQUS  built-in  von  Mises  and  Hill  material  models  for  the  matrix. 
As  an  example,  in  Fig.  10  are  shown  the  results  of  these  FE  calcula¬ 
tions  in  terms  of  the  von  Mises  effective  stress  Ze  versus  the  von 
Mises  effective  strain  Ee  for  each  material.  In  the  calculations  illus¬ 
trated  in  Fig.  10,  the  displacements  were  prescribed  such  that  a 
constant  value  of  the  strain  triaxiality  (TE  =  1.5)  was  maintained 
and  the  major  stress  was  oriented  along  the  e3  direction  (Z33  > 
Zn  >  0).  Both  materials  B  and  C  have  the  same  yield  strength  along 
the  in-plane  direction  ei.  However,  since  the  matrix  in  material  C  is 
softest  in  the  e3  direction,  there  is  more  plastic  flow  and,  hence, 
more  void  growth  than  in  material  B,  for  which  the  e3  direction 
is  the  hardest.  Consequently,  material  C  yields  first  at  at  a  lower 
stress  level  and  its  overall  softening  is  much  more  pronounced 
than  in  material  B.  Comparison  between  the  isotropic  stress-strain 
curve  and  the  anisotropic  ones  illustrate  very  clearly  the  effect  of 
the  directionality  of  the  plastic  flow  of  the  matrix  on  yielding  of 
the  void  matrix  aggregate.  Overall  softening  in  the  isotropic  mate¬ 
rial  is  less  pronounced  than  in  material  C  since,  while  all  three 
materials  have  the  same  yield  strength  in  the  ei  direction,  material 
A  has  equal  yield  strengths  in  the  ei  and  e3  directions  whereas 
material  C  has  a  lower  yield  strength  in  the  e3  direction. 

Figs.  11-13  show  the  finite  element  results  versus  Eq.  (86) 
where  three  different  materials  (corresponding  to  material  types 
A,  B  and  C  with  k  =  0.3098  such  that  the  yield  strengths  in  tension 
are  greater  than  in  compression)  are  each  used  as  the  matrix  mate¬ 
rial  in  a  unit  cell  containing  a  spherical  void  (1%  void  volume  frac¬ 
tion).  The  figures  show  the  normalized  macroscopic  von  Mises 
effective  stress  Ze/a\  versus  the  normalized  macroscopic  mean 
stress  Zm/(j[;  the  numerical  values  of  the  anisotropy  coefficients 
are  given  in  Table  1.  The  macroscopic  yield  points  for  the  finite  ele¬ 
ment  calculations  were  taken  to  correspond  with  the  points  of 
maximum  macroscopic  effective  stress.  Results  corresponding  to 
both  positive  and  negative  values  of  the  third  invariant  of  the 
stress  deviator  fl  are  reported  in  the  figures. 

Fig.  1 1  shows  a  comparison  between  the  finite  element  results 
versus  the  theoretical  yield  curves  given  by  Eq.  (86)  when  the 
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Fig.  10.  Typical  von  Mises  effective  stress  versus  von  Mises  effective  strain  curves  corresponding  to  FE  calculations  at  TE=  1.5  and  with  the  major  stress  in  the  e3  direction.  For 
all  materials  the  yield  stress  in  tension  and  compression  are  equal  (/<  =  0)  with  material  type  A  being  isotropic  (< t\,c  =  crj0),  material  type  B  having  a  through-thickness  yield 
strength  greater  than  the  in-plane  yield  strength  (er  V  <  v?)  and  material  type  C  having  a  through-thickness  yield  strength  lower  than  the  in-plane  yield  strength 
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Fig.  11.  Analytical  yield  curves  according  to  the  Benzerga  and  Besson  (2001 )  criterion  and  axisymmetric  FE  results  corresponding  to  an  anisotropic  matrix  (i.e.,  material  type 
B  in  Table  1  with  k  =  0).  The  void  volume  fraction  is/0  =  0.01. 


matrix  materials  are  anisotropic  (the  anisotropy  coefficients  corre¬ 
spond  to  material  type  B  of  Table  1)  and  does  not  display  tension 
compression  asymmetry  (/<  =  0).  Thus,  Eq.  (86)  reduces  to  the  ana¬ 
lytical  criterion  of  Benzerga  and  Besson  (2001).  As  in  the  case  of 
Gurson  (1977),  the  Benzerga  and  Besson  (2001 )  yield  criterion  does 
not  account  for  the  effect  of  the  third  invariant  of  the  stress  devi¬ 
ator;  yet,  the  agreement  with  FE  results  is  still  quite  good. 

Figs.  12  and  13  show  FE  results  together  with  the  analytical 
yield  curves  corresponding  to  materials  having  the  yield  strength 
in  the  through-thickness  direction  greater  than  the  yield  strength 
in  any  in-plane  direction  and  vice  versa,  respectively  (i.e.,  materials 
type  B  and  C,  respectively,  from  Table  1 ).  To  examine  the  combined 
effects  of  anisotropy  and  tension-compression  asymmetry,  axi¬ 
symmetric  calculations  were  performed  with  the  through-thick¬ 
ness  direction,  e3,  corresponding  to  either  the  major  stress 


C/J  >  0)  or  to  the  minor  stress  (J 3  <  0).  The  agreement  between 
the  analytical  yield  curves  according  to  the  developed  criterion 
(Eq.  (86))  and  the  finite  element  results  obtained  from  the  unit  cell 
calculations  is  satisfactory. 

Finally,  Fig.  14  shows  a  comparison  between  the  FE  results  and 
the  anisotropic  yield  criterion  of  Eq.  (86)  for  three  different  mate¬ 
rials,  each  having  the  through-thickness  direction  e3  the  hardest 
but  displaying  various  levels  of  tension-compression  asymmetry 
(/<  =  0,  k  >  0  and  k  <  0).  The  FE  data  corresponds  to  axisymmetric 
loading  with  the  minor  stress  being  in  the  through-thickness  direc¬ 
tion  (i.e.,  J3  <  0).  The  model  accurately  predicts  that  for  the  same 
yield  strength  ratio,  the  softest  material  corresponds  to  k  >  0  while 
the  hardest  corresponds  to  k  <  0  (see  also  Eqs.  (70),  (71)  and  (77) 
which  give  the  expressions  of  the  intercepts  of  the  yield  locus  with 
the  hydrostatic  and  deviatoric  axes  as  functions  of  the  strength- 
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Fig.  12.  Analytical  yield  curves  according  to  the  proposed  criterion  (Eq.  (86))  and  axisymmetric  FE  results  for  both  jf  >  0  and  il<  0.  The  anisotropic  matrix  material  is 
hardest  in  the  through-thickness  direction  (i.e.,  material  type  B  in  Table  1)  and  displays  tension-compression  asymmetry  with  the  tensile  yield  strengths  larger  than  the 
compressive  ones  ( k  =  0.3098).  The  void  volume  fraction  is /0  =  0.01. 
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Fig.  13.  Analytical  yield  curves  according  to  the  proposed  criterion  (Eq.  (86))  and  axisymmetric  FE  results  for  both jf  >  0  and /f  <  0.  The  anisotropic  matrix  material  is  softest 
in  the  through-thickness  direction  (i.e.,  material  type  C  in  Table  1)  and  displays  tension-compression  asymmetry  with  the  tensile  yield  strengths  larger  than  the  compressive 
ones  (/<  =  0.3098).  The  void  volume  fraction  is/0  =  0.01. 


differential  parameter  k).  Once  again,  the  comparison  between  fi¬ 
nite  element  data  and  the  yield  curves  shows  quite  good 
agreement. 


7.  Concluding  remarks 

Yielding  of  porous  materials  in  which  the  matrix  is  incompress¬ 
ible,  anisotropic,  and  displays  tension-compression  asymmetry 
has  been  studied.  An  analytical  yield  criterion  has  been  developed 
by  extending  Gurson’s  (1977)  analysis  of  the  hollow  sphere  to  the 
case  when  the  matrix  plastic  behavior  is  described  by  the  aniso¬ 
tropic  yield  criterion  of  Cazacu  et  al.  (2006).  The  classical  velocity 
field  of  Rice  and  Tracey  (1969)  has  been  used  in  the  upper  bound 
calculations.  This  velocity  field  is  the  simplest  among  the  fields 
that  meet  the  requirements  of  the  Hill-Mandel  lemma,  i.e.,  incom¬ 


pressibility  and  compatibility  with  uniform  strain  rate  boundary 
conditions.  There  are  several  studies  that  motivate  the  examina¬ 
tion  of  more  complex  microscopic  flow  fields  (e.g.,  Needleman, 
1972;  Koplik  and  Needleman,  1988).  Specifically,  it  is  shown  that 
part  of  the  matrix  might  not  attain  plastic  yield  and  that  elastic 
unloading  may  occur  in  certain  subdomains.  Nevertheless,  even 
using  the  simplest  velocity  field,  fresh  difficulties  were  encoun¬ 
tered  when  estimating  the  local  plastic  dissipation  associated  with 
the  Cazacu  et  al.  (2006)  yield  criterion.  These  difficulties  are  related 
to  the  tension-compression  asymmetry  and  anisotropy  of  the  ma¬ 
trix  response.  Thus,  the  plastic  multiplier  rate  has  multiple 
branches  (see  Eq.  (28))  and  depends  on  each  of  the  principal  values 
of  the  local  rate  of  deformation,  d. 

Certain  approximations  were  introduced  in  order  to  obtain  the 
analytical,  closed-form  expression  of  the  overall  plastic  potential 
(Eq.  (68)).  The  derived  criterion  is  anisotropic,  exhibits  tension- 
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Fig.  14.  Anisotropic,  axisymmetric  finite  element  yield  points  versus  analytic  curves  when  the  initial  void  volume  fraction  is/0  =  0.01  and  the  matrix  material  is  hardest  in  the 
through-thickness  direction  (i.e.,  material  type  B  in  Table  1).  The  curves  shown  are  for  k  =  0  (no  strength  differential),  k  =  0.3098  (yield  strengths  in  tension  greater  than  in 
compression)  and  k  =  -0.3098  (yield  strengths  in  tension  less  than  in  compression)  with  Z33  < 


compression  asymmetry  due  to  the  characteristics  of  the  plastic 
flow  of  the  matrix  and  is  pressure-sensitive  due  to  the  presence 
of  voids.  It  is  worth  noting  that  this  criterion  depends  on  all  stress 
invariants  as  well  as  the  mixed  invariants  between  stress  and 
structural  tensors.  Comparison  between  the  theoretical  predictions 
using  this  criterion  and  results  of  finite  element  cell  calculations 
show  an  overall  good  agreement. 

Although  the  expression  of  the  proposed  criterion  for  the  por¬ 
ous  material  is  similar  to  that  of  Benzerga  and  Besson  (2001 ),  there 
are  distinct  differences. 

•  The  Hill  (1948)  equivalent  stress  is  replaced  by 

Ie  =  rhyjYd=i(\Zi\  ~  k^')2>  the  equivalent  stress  associated  to 
the  CPB06  anisotropic  criterion,  which  depends  on  all  principal 
values  of  the  stress  deviator  £',  the  anisotropy  coefficients  Ly 
and  on  the  ratio  between  the  uniaxial  yield  in  tension  g\  and 
the  uniaxial  yield  in  compression  of  of  the  matrix  through  the 
parameter  m  (see  Eq.  (16)  for  the  definition  of  the  constant  m). 

•  It  involves  a  new  coefficient  h  such  that  for  tensile  hydrostatic 
loading,  yielding  of  the  void-matrix  aggregate  occurs  when 

Zm  =  -a\  ln(f )  while  for  compressive  hydro¬ 

static  loading  yielding  is  at  Zm  =  <r\ ln(/)-  Thus, 
for  arbitrary  loadings  the  effective  yield  locus  is  no  longer  sym¬ 
metric  with  respect  to  the  vertical  axis  £m  =  0,  as  it  is  in  the  case 
of  Gurson  (1977)  or  Benzerga  and  Besson’s  (2001)  yield  locus. 

If  there  is  no  difference  in  response  between  the  yield  in  tension 
and  compression,  Ie  reduces  to  Hill’s  equivalent  stress  and  the 
coefficient  h  reduces  to  the  expression 

h  =  8/15)6  +h  +  h)  +  (4/5)(/4  +  /5  +  40,  (87) 

which  coincides  with  the  hydrostatic  factor  in  Benzerga  and  Bes¬ 
son’s  (2001)  criterion  (since  for  k  =  0,  the  quadratic  form  of  Cazacu 
et  al.’s  (2006)  anisotropic  criterion  reduces  to  Hill’s  (1948)  aniso¬ 
tropic  criterion);  hence,  the  proposed  criterion  (Eq.  (68))  reduces 
to  the  criterion  of  Benzerga  and  Besson  (2001).  In  the  absence  of 
voids,  the  proposed  criterion  reduces  to  Cazacu  et  al.’s  (2006)  aniso¬ 
tropic  yield  criterion.  The  accuracy  of  the  analytical  criterion  was 


assessed  through  comparison  with  finite-element  cell  calculations. 
To  improve  the  agreement,  the  proposed  analytic  yield  criterion 
(Eq.  (68))  was  modified  to  include  additional  parameters,  qi}  as  done 
by  Tvergaard  (1981),  Tvergaard  and  Needleman  (1984)  in  the  case 
of  the  Gurson  (1977)  yield  criterion.  In  this  manner,  for  k  =  0  and 
L  =  1  (a  von  Mises  matrix),  the  criterion  reduces  to  the  GTN  model, 
while  for  k  non-zero  and  L  isotropic  it  reduces  to  Cazacu  and  Stew¬ 
art  (2009).  The  agreement  between  the  theoretical  predictions 
using  this  criterion  (Eq.  (86))  and  the  results  of  finite  element  cell 
calculations  are  quite  good. 
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Appendix  A.  Calculation  of  the  hydrostatic  parameter 

The  hydrostatic  parameter,  h,  is  given  as 

h2  =  h(d  :  L  :  d)s(r)  =  —  J  J  (d  :  L  :  dj  sin 0d0d0,  (88) 

where 

d  =  -2er  (g)  er  +  ee  <8>  e6  +  e^  ®  e0.  (89) 

The  basis  vectors  for  the  spherical  coordinate  system  (er,  e0,  e^)  are 
expressed  in  terms  of  the  basis  vectors  for  the  Cartesian  coordinate 
system  associated  with  the  material  axes  of  symmetry  (ei,  e2,  e3)  as 
follows: 

er  =  (sin  0  cos  0)e  i  +  (sin  0  sin  0)e2  +  cos  0e3 , 
e0  =  (cos  6  cos  0)e i  +  (cos  6  sin  0)e2  -  sin  0e3l  (90) 

=  -sin^ei  +cos  0e2, 

where  0  e  [0,7i]  and  </>  e  [0,27i]. 

In  order  to  evaluate  the  integral  expression  in  Eq.  (88),  the  ten¬ 
sor  d  will  be  transformed  from  the  spherical  coordinate  system  to 
the  material  coordinate  system.  Thus, 


J.B.  Stewart,  0.  Cazacu/ International  Journal  of  Solids  and  Structures  48  (20U)  357-373 


373 


d(i,2,3)  —  Rd(r)0j(£)Rr, 


where 

"cos^sin  9 

cos  (j)  cos  9 

-  sin  0 

R  = 

sin  (p  sin  9 

sin  0  cos  0 

COS(/> 

cos9 

-sine 

0 

such  that 

(91) 


(92) 


d 


(1,2,3) 


1  -3sin20cos2</>  -3  sin2  0sin  (p  cos  </> 

-3  sin2  9  sin  </>  cos  </>  1-3  sin2  9  sin2  </> 

-3  sin  9  cos  9  cos  </>  -3sin0cos0sin</> 


-3sin0cos0cos</> 
-3sin0cos0sin</>  • 
1-3  cos2  9 

(93) 


Note  that 


a  :  L  :  a  =  /id,,  +  hdj2  +  hdj3  +  2 (l4d223  +  h~d23  +  i6d?2) .  (94) 


Therefore,  evaluating  the  integral  expression  in  Eq.  (88)  reduces  to 
evaluating  the  surface  integral  of  the  squared  components  of  d(ij2>3). 
Referring  to  Eq.  (93)  yields 

(^ll)s(r)  —  (^22 )s(r)  =  (^33 )s(r)  =  r  > 

3  (95) 

(^23  )s(r)  =  (^13  )s(r)  =  (^12  )s(r)  =  5  • 


Finally,  an  expression  for  the  hydrostatic  parameter,  h ,  is  arrived  at 
by  combining  Eqs.  (88),  (94)  and  (95).  Thus, 


hz=h 


^  (^1  +^2+^3^  +  ^  ^4  +  h  + 


(96) 
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